Common solutions for creating maps usually involves GIS software, such as ArcGIS, QGIS, eSpatial, etc., which allow to visually prepare a map, in the same approach as one would prepare a poster or a document layout. On the other hand, R has also advanced spatial capabilities and can be used to draw maps. As we mentioned earlier for the graphics, R allows you to create map highly customization which makes it a solution increasingly attractive.
Several packahges will be used today. It is generally a good practice to group the package you will use at the beginning of your script:
library(maptools)
library(rgdal)
library (sp)
library(raster)
library(ggplot2)
library(sf)
library(rnaturalearth)
library(rnaturalearthdata)
library(ggspatial)
library(rgbif)
library(mapr)
library(marmap)
library(leaflet)
#Sys.setlocale("LC_TIME", "English")
#Sys.setlocale(category = "LC_ALL", "Chinese (Traditional)_Taiwan.950")
Getting a simple map is pretty straightforward in R using the
packagemaptools
data(wrld_simpl)
# NOT RUN: ? wrld_simpl
plot(wrld_simpl)
Now we know the options to modify a plot, therefore you can customize this map with some elements we mentioned earlier:
Try to zoom on Taiwan: you will quickly realize that we have a pretty
bad resolution a finer scale, see ?wrld_simpl
gpxIf you like running, hiking, walking, or being lazy in a park - many
devices now can export gpx format. You an read the file
directly as lines (i.e, tracks), points (i.e,
track_points), and a few other formats you can find in the
help for readOGR from the famous package
rgdal: bindings for the ‘Geospatial’ Data Abstraction
Library. Note that the same could be accomplished using the
sf package and st_readfunction.
OGR data source with driver: GPX
Source: "C:\Users\Vdenis\OneDrive\Documents\Class\OCEAN50998\Website\OCEAN5098\8a58a0ea3eb13bdbea53dfe957e0e06a7c315f2c\Data\run.gpx", layer: "tracks"
with 1 features
It has 12 fields
plot(run1, main='Line') # my running activity
run2 <- readOGR(dsn="Data/run.gpx",layer="track_points")
OGR data source with driver: GPX
Source: "C:\Users\Vdenis\OneDrive\Documents\Class\OCEAN50998\Website\OCEAN5098\8a58a0ea3eb13bdbea53dfe957e0e06a7c315f2c\Data\run.gpx", layer: "track_points"
with 931 features
It has 26 fields
plot(run2, main='Points')
dev.off()
RStudioGD
2
gpxThe reverse can be accomplish with the function writeOGR
to export objects as SpatialPointsDataFrame,
SpatialLinesDataFrame, or
SpatialPolygonsDataFrame (as defined in the
sp package). This will write an ArcGIS compatible
shapefile, but many different formats are available by
specifying the correct driver.Now that we know how to get a basic map in
R, let’s look at how we can export and import data.
writeOGR(wrld_simpl, dsn="Data", layer = "world_test", driver = "ESRI Shapefile", overwrite_layer = TRUE)
shp filesNow we could open world_test.shp in ArcGIS (or others),
but we can also import shapefile (.shp) back into R, let’s use that same
file
world_shp <- readOGR(dsn = "Data",layer = "world_test")
OGR data source with driver: ESRI Shapefile
Source: "C:\Users\Vdenis\OneDrive\Documents\Class\OCEAN50998\Website\OCEAN5098\8a58a0ea3eb13bdbea53dfe957e0e06a7c315f2c\Data", layer: "world_test"
with 246 features
It has 11 fields
plot(world_shp)

Creating spatial data from scratch in R seems a little bit confusing at first, but once you got the logic behind, it get easier. As we learned before for graphics you can add vector based object (points, Lines, polygones) to your map. A transformation of the numerical x and y values applied in this case to recognize those values as spatial coordinate.
SpatialPointsDataFrame for plotting pointsplot(wrld_simpl,xlim=c(115,128) ,ylim=c(19.5,27.5),col='#D2B48C',bg='lightblue') # TW map
coords <- matrix(c(121.537290,120.265541, 25.021335, 22.626524),ncol=2) # NTU and SYS univ.
coords <- coordinates(coords) # assign values as spatial coordinates
spoints <- SpatialPoints(coords) # create SpatialPoints
df <- data.frame(location=c("NTU","SYS")) # create a dataframe
spointsdf <- SpatialPointsDataFrame(spoints,df) # create a SpatialPointsDataFrame
plot(spointsdf,add=T,col=c('black','black'),pch=19,cex=2.2) # plot it on our map
text(121,24, 'TAIWAN', cex=1)

SpatialLinesDataFrame for plotting lines: let’s travel
in Canada in the province of Saskatchewanplot(wrld_simpl,xlim=c(-130,-60),ylim=c(45,80),col='#D2B48C',bg='lightblue')
coords <- matrix(c(-110,-102,-102,-110,-110,60,60,49,49,60),ncol=2)
l <- Line(coords)
ls <- Lines(list(l),ID="1")
sls <- SpatialLines(list(ls))
df <- data.frame(province="Saskatchewan")
sldf <- SpatialLinesDataFrame(sls,df)
plot(sldf,add=T,col='#3d2402', cex=2)
text(-114, 55, 'Saskatchewan', srt=90, cex=0.7)
text(-114, 63, 'CANADA', cex=1)

SpatialPolygonesDataFrame for plotting polygonesplot(wrld_simpl,xlim=c(-130,-60),ylim=c(45,80),col='#D2B48C',bg='lightblue')
coords <- matrix(c(-110,-102,-102,-110,-110,60,60,49,49,60),ncol=2)
p <- Polygon(coords)
ps <- Polygons(list(p),ID="1")
sps <- SpatialPolygons(list(ps))
df <- data.frame(province="Saskatchewan")
spdf <- SpatialPolygonsDataFrame(sps,df)
plot(spdf,add=T,col='#45220d')
text(-114, 55, 'Saskatchewan', srt=90, cex=0.7)
text(-114, 63, 'CANADA', cex=1)
text(-103, 46, 'UNITED STATES', cex=1)
text(-40, 78, 'GREENLAND', cex=1)
text(-35, 55, 'Atlantic Ocean', cex=1, col='#071238')

raster fct.Spatial data: download, unzip, and import spatial shape (don’t forget to set up your working directory). Those spatial data are high resolution spatial data at the country level.
Using the rasterpackage and the function
getData, you can easily download directly polygones
(vectorized) shape from GADM (Global Administrative
Areas, or other sources). You will simply need the country code such as
TWN for Taiwan, JPN for Japan, etc. In
addition, the argument level indicates the level you are
trying to get. For provinces (first level subdivision))
level=1, for the overall country level
level=0. Check ?getDta.
TWN1 <- getData('GADM', country="TWN", level=0) # data Taiwan
JPN <- getData('GADM', country="JPN", level=0) # data Japan
class(TWN1) # those datasets are SpatialPolygonsDataFrame
[1] "SpatialPolygonsDataFrame"
attr(,"package")
[1] "sp"
par(mfrow = c(1, 2))
plot(TWN1,axes=T,bg=colors()[431],col='grey')
plot(JPN,axes=T,bg=colors()[431],col='grey')

Don’t forget to close your graphical window:
dev.off()
RStudioGD
2
Simply set up xlim and ylim
TWN2 <- getData('GADM', country="TWN", level=1)
TWN2$NAME_1
[1] "Fujian" "Kaohsiung" "New Taipei" "Taichung" "Tainan"
[6] "Taipei" "Taiwan"
You can see the list of counties is not complete. It is often the case when using GADM that some ‘details’ will be missing. Let’s visualize one county for which we have the data shape.
plot(TWN1,col="grey",xlim=c(119,122.5), ylim=c(21.5,25.5), bg=colors()[431], axes=T)
KAO <- TWN2[TWN2$NAME_1=="Kaohsiung",]
plot(KAO,col="grey 33",add=TRUE)

Note the use of add=TRUE in the latest plot
function.
# base map
plot(TWN1,col="grey",xlim=c(121,122), ylim=c(24,25.5), bg=colors()[431], axes=T)
# adding spatial polygones
TAI <- TWN2[TWN2$NAME_1=="Taipei" | TWN2$NAME_1=="New Taipei" ,]
plot(TAI,col="black",add=TRUE)
# adding spatial points
coords <- matrix(cbind(lon=c(121.2,121.55,121.8),lat=c(25,25.19,24.5)),ncol=2)
coords <- coordinates(coords)
spoints <- SpatialPoints(coords)
df <- data.frame(location=c("City 1","City 2","City 3"),pop=c(138644,390095,34562))
spointsdf <- SpatialPointsDataFrame(spoints,df)
scalefactor <- sqrt(spointsdf$pop)/sqrt(max(spointsdf$pop))
plot(spointsdf,add=TRUE,col='white',pch=1,cex=scalefactor*3,lwd=2)
# adding a location of NTU (not spatial point here)
points(121.537290,25.021335, type="p", pch=18, col='white', cex=1.5)
# adding text
text(121.53,24.921335,"NTU", col='white', font=2)
# adding scale
maps::map.scale(x=120, y=25.4)
# adding north arrow
GISTools::north.arrow(xb=120.3,yb=24.7, len=0.06, lab='N')

Note the use of the synthax maps::map.scale. It means
you wanna use the function map.scale from the package
maps. It avoid passing by library(maps), thus
avoiding possible conflict among functions.
Each country will usually have open data platform with country shape data. Those data are usually more accurate and should be preferred. In Taiwan, shapefile of counties are available for download here. You could get even get a finer resolution at the township level if you wish (large file).
First prepare your system to receive traditional Chinese characters:
Sys.setlocale(category = "LC_ALL", "Chinese (Traditional)_Taiwan.950")
[1] "LC_COLLATE=Chinese (Traditional)_Taiwan.950;LC_CTYPE=Chinese (Traditional)_Taiwan.950;LC_MONETARY=Chinese (Traditional)_Taiwan.950;LC_NUMERIC=C;LC_TIME=Chinese (Traditional)_Taiwan.950"
R version 4.2.1 (2022-06-23 ucrt)
Platform: x86_64-w64-mingw32/x64 (64-bit)
Running under: Windows 10 x64 (build 19043)
Matrix products: default
locale:
[1] LC_COLLATE=Chinese (Traditional)_Taiwan.950
[2] LC_CTYPE=Chinese (Traditional)_Taiwan.950
[3] LC_MONETARY=Chinese (Traditional)_Taiwan.950
[4] LC_NUMERIC=C
[5] LC_TIME=Chinese (Traditional)_Taiwan.950
system code page: 65001
attached base packages:
[1] stats graphics grDevices utils datasets methods
[7] base
other attached packages:
[1] psych_2.2.5 MVPARTwrap_0.1-9.2
[3] mvpart_1.6-2 mvabund_4.2.1
[5] devtools_2.4.4 usethis_2.1.6
[7] caret_6.0-93 randomForest_4.7-1.1
[9] rattle_5.5.1 bitops_1.0-7
[11] tibble_3.1.7 rpart_4.1.16
[13] tree_1.0-42 e1071_1.7-11
[15] gridExtra_2.3 fpc_2.2-9
[17] factoextra_1.0.7 qgraph_1.9.2
[19] gclus_1.3.2 cluster_2.1.3
[21] ade4_1.7-19 ape_5.6-2
[23] vegan_2.6-2 permute_0.9-7
[25] marmap_1.0.6 mapr_0.5.2
[27] rgbif_3.7.2 ggspatial_1.1.6
[29] rnaturalearthdata_0.1.0 rnaturalearth_0.1.0
[31] sf_1.0-8 raster_3.5-29
[33] rgdal_1.5-32 maptools_1.1-4
[35] sp_1.5-0 gganimate_1.0.7
[37] animation_2.7 nlme_3.1-157
[39] lmerTest_3.1-3 lme4_1.1-30
[41] readr_2.1.2 yarrr_0.1.5
[43] circlize_0.4.15 BayesFactor_0.9.12-4.4
[45] Matrix_1.4-1 coda_0.19-4
[47] jpeg_0.1-9 interactions_1.1.5
[49] car_3.1-0 carData_3.0-5
[51] MASS_7.3-57 corrplot_0.92
[53] Hmisc_4.7-1 Formula_1.2-4
[55] survival_3.3-1 leaflet_2.1.1
[57] forcats_0.5.1 hexbin_1.28.2
[59] stringr_1.4.0 tidyr_1.2.0
[61] ggplot2_3.3.6 lattice_0.20-45
[63] readxl_1.4.0 dplyr_1.0.9
[65] scales_1.2.0
loaded via a namespace (and not attached):
[1] corpcor_1.6.10 class_7.3-20
[3] ps_1.7.1 foreach_1.5.2
[5] rprojroot_2.0.3 crayon_1.5.1
[7] fBasics_4021.92 backports_1.4.1
[9] rlang_1.0.3 nloptr_2.0.3
[11] callr_3.7.0 bit64_4.0.5
[13] glue_1.6.2 parallel_4.2.1
[15] processx_3.6.1 classInt_0.4-7
[17] tidyselect_1.1.2 ggpubr_0.4.0
[19] xtable_1.8-4 MatrixModels_0.5-0
[21] magrittr_2.0.3 evaluate_0.16
[23] ncdf4_1.19 cli_3.3.0
[25] rstudioapi_0.13 miniUI_0.1.1.1
[27] whisker_0.4 bslib_0.3.1
[29] jtools_2.2.0 maps_3.4.0
[31] GISTools_0.7-4 shiny_1.7.2
[33] xfun_0.31 pkgbuild_1.3.1
[35] ggrepel_0.9.1 listenv_0.8.0
[37] installr_0.23.2 png_0.1-7
[39] future_1.27.0 ipred_0.9-13
[41] withr_2.5.0 tweedie_2.3.5
[43] plyr_1.8.7 cellranger_1.1.0
[45] hardhat_1.2.0 pROC_1.18.0
[47] pillar_1.8.1 GlobalOptions_0.1.2
[49] cachem_1.0.6 fs_1.5.2
[51] flexmix_2.3-18 kernlab_0.9-31
[53] vctrs_0.4.1 pbivnorm_0.6.0
[55] ellipsis_0.3.2 oai_0.3.2
[57] generics_0.1.3 lava_1.6.10
[59] urltools_1.7.3 tools_4.2.1
[61] foreign_0.8-82 munsell_0.5.0
[63] tweenr_2.0.0 proxy_0.4-27
[65] fastmap_1.1.0 compiler_4.2.1
[67] pkgload_1.3.0 abind_1.4-5
[69] httpuv_1.6.5 sessioninfo_1.2.2
[71] rgeos_0.5-9 prodlim_2019.11.13
[73] rpart.plot_3.1.1 deldir_1.0-6
[75] utf8_1.2.2 later_1.3.0
[77] recipes_1.0.1 wk_0.6.0
[79] jsonlite_1.8.0 pbapply_1.5-0
[81] lazyeval_0.2.2 promises_1.2.0.1
[83] latticeExtra_0.6-30 checkmate_2.1.0
[85] rmarkdown_2.14 textshaping_0.3.6
[87] statmod_1.4.37 Rtsne_0.16
[89] pander_0.6.5 igraph_1.3.4
[91] numDeriv_2016.8-1.1 yaml_2.3.5
[93] systemfonts_1.0.4 prabclus_2.3-2
[95] conditionz_0.1.0 htmltools_0.5.2
[97] memoise_2.0.1 lavaan_0.6-12
[99] profvis_0.3.7 modeltools_0.2-23
[101] viridisLite_0.4.0 digest_0.6.29
[103] assertthat_0.2.1 mime_0.12
[105] adehabitatMA_0.3.14 units_0.8-0
[107] RSQLite_2.2.16 future.apply_1.9.0
[109] remotes_2.4.2 data.table_1.14.2
[111] urlchecker_1.0.1 blob_1.2.3
[113] ragg_1.2.2 ggsci_2.9
[115] splines_4.2.1 labeling_0.4.2
[117] broom_1.0.0 hms_1.1.1
[119] colorspace_2.0-3 base64enc_0.1-3
[121] mnormt_2.1.0 shape_1.4.6
[123] nnet_7.3-17 sass_0.4.1
[125] Rcpp_1.0.9 mclust_5.4.10
[127] bookdown_0.28 mvtnorm_1.1-3
[129] fansi_1.0.3 tzdb_0.3.0
[131] parallelly_1.32.1 ModelMetrics_1.2.2.2
[133] R6_2.5.1 grid_4.2.1
[135] crul_1.2.0 lifecycle_1.0.1
[137] formatR_1.12 zip_2.2.0
[139] ggsignif_0.6.3 curl_4.3.2
[141] minqa_1.2.4 jquerylib_0.1.4
[143] robustbase_0.95-0 distill_1.4
[145] glasso_1.11 RColorBrewer_1.1-3
[147] iterators_1.0.14 gower_1.0.0
[149] htmlwidgets_1.5.4 triebeard_0.3.0
[151] purrr_0.3.4 crosstalk_1.2.0
[153] terra_1.6-7 labdsv_2.0-1
[155] timeSeries_4021.104 mgcv_1.8-40
[157] globals_0.16.0 htmlTable_2.4.1
[159] codetools_0.2-18 lubridate_1.8.0
[161] gtools_3.9.3 leaflet.providers_1.9.0
[163] prettyunits_1.1.1 gtable_0.3.0
[165] DBI_1.1.3 stats4_4.2.1
[167] httr_1.4.3 highr_0.9
[169] KernSmooth_2.23-20 stringi_1.7.6
[171] progress_1.2.2 spatial_7.3-15
[173] reshape2_1.4.4 farver_2.1.1
[175] uuid_1.1-0 diptest_0.76-0
[177] fdrtool_1.2.17 magick_2.7.3
[179] timeDate_4021.104 xml2_1.3.3
[181] boot_1.3-28 s2_1.1.0
[183] interp_1.1-3 DEoptimR_1.0-11
[185] bit_4.0.4 pkgconfig_2.0.3
[187] rstatix_0.7.0 knitr_1.39
[189] downlit_0.4.2 httpcode_0.3.0
You can go to the website and downaload data. Alternatively, using
the code source of the page your can target directly the the
url of teh dataset you want to download. This is a zip file
that you want to temporary store and extract in a specific directory
url <- 'https://data.moi.gov.tw/MoiOD/System/DownloadFile.aspx?DATA=72874C55-884D-4CEA-B7D6-F60B0BE85AB0'
path1 <- tempfile(fileext = ".zip")
if (file.exists(path1)) 'file alredy exists' else download.file(url, path1, mode="wb")
zip::unzip(zipfile = path1,exdir = 'Data')
Again, the function readOGR reads the .shp file. Other
arguments are for the file encoding and to deal with the Chinese
characters because my system is in English.
taiwan <- readOGR('Data/COUNTY_MOI_1090820.shp', use_iconv=TRUE, encoding='UTF-8')
OGR data source with driver: ESRI Shapefile
Source: "C:\Users\Vdenis\OneDrive\Documents\Class\OCEAN50998\Website\OCEAN5098\8a58a0ea3eb13bdbea53dfe957e0e06a7c315f2c\Data\COUNTY_MOI_1090820.shp", layer: "COUNTY_MOI_1090820"
with 22 features
It has 4 fields
Make the plot:
This is the official map of TAIWAN, including all remote territories such as Taiping and Diaoyutai islands. County names are here accessible using:
taiwan$COUNTYNAME
[1] "連江縣" "宜蘭縣" "彰化縣" "南投縣" "雲林縣" "基隆市" "臺北市"
[8] "新北市" "臺中市" "臺南市" "桃園市" "苗栗縣" "嘉義市" "嘉義縣"
[15] "金門縣" "高雄市" "臺東縣" "花蓮縣" "澎湖縣" "新竹市" "新竹縣"
[22] "屏東縣"
taiwan$COUNTYENG
[1] "Lienchiang County" "Yilan County" "Changhua County"
[4] "Nantou County" "Yunlin County" "Keelung City"
[7] "Taipei City" "New Taipei City" "Taichung City"
[10] "Tainan City" "Taoyuan City" "Miaoli County"
[13] "Chiayi City" "Chiayi County" "Kinmen County"
[16] "Kaohsiung City" "Taitung County" "Hualien County"
[19] "Penghu County" "Hsinchu City" "Hsinchu County"
[22] "Pingtung County"
ggplot2mapsAs we mentioned earlier, the package ggplot2 implements
the grammar of graphics in R. While ggplot2 is becoming the
standard for R graphs, it does not handle spatial data specifically. The
current state-of-the-art of spatial objects in R relies on Spatial
classes defined in the package sp, but the new package
sf has implemented the “simple feature” standard, and is
steadily taking over sp. Recently, the package
ggplot2 has allowed the use of simple features from the
package sf as layers in a graph (since the version 3.0.0 of
ggplot2). The combination of ggplot2 and
sf therefore enables to programmatically create maps, using
the grammar of graphics, just as informative or visually appealing as
traditional GIS software.
Note that the conversion from sp to sf is
acheivable using the function st_as_sf from the package
sf.
# not run
# wrld_simpl <- st_as_sf(wrld_simpl)
After loading the basic packages necessary for all maps,
i.e. ggplot2 and sf. We can change the theme
of our plot. The classic dark-on-light theme for ggplot2
(theme_bw) is appropriate for maps (default:
theme_grey()):
The package rnaturalearth provides a map of countries of
the entire world. Use ne_countries to pull country data and
choose the scale (rnaturalearthhires is necessary for
scale = "large"). The function can return sp
classes (default) or directly sf classes, as defined in the
argument returnclass:
world <- ne_countries(scale = "medium", returnclass = "sf")
class(world)
[1] "sf" "data.frame"
In addition to country shape polygones, this dataset contains information on the population in every country.
ggplot mapFirst, let us start with creating a base map of the world using
ggplot2. This base map will then be extended with different
map elements, as well as zoomed in to an area of interest. We can check
that the world map was properly retrieved and converted into an
sf object, and plot it with ggplot2. This time
we will indicate that the geometry of our plot is constrained by a
sf object:
In your map of the world, the plot panel is expanded beyond the size of the earth (you can see that the graticule lines end before the edge of the plot panel), and hence no axis ticks are drawn. One way to solve the issue is to turn off the expansion while defining the coordinates.
As before, the layers are added one at a time in a
ggplot call, so the order of each layer is very important.
All data will have to be in an sf format to be used by
ggplot2; data in other formats (e.g. classes from
sp) will be manually converted to sf classes
if necessary.
ggtitle,
xlab, ylabA title and a subtitle can be added to the map using the function
ggtitle, passing any valid character string (e.g. with
quotation marks) as arguments. Axis names are absent by default on a
map, but can be changed to something more suitable (e.g. “Longitude” and
“Latitude”), depending on the map:
ggplot(data = world) +
geom_sf() +
coord_sf(expand = FALSE) +
xlab("Longitude") + ylab("Latitude") +
ggtitle("World map", subtitle = paste0("(", length(unique(world$name)), " countries)"))

geom_sfIn many ways, sf geometries are no different than
regular geometries, and can be displayed with the same level of control
on their attributes. Here is an example with the polygons of the
countries filled with a green color (argument fill), using
black for the outline of the countries (argument
color):
The package ggplot2 allows the use of more complex color
schemes, such as a gradient on one variable of the data. Here is another
example that shows the population of each country. In this example, we
use the “viridis” colorblind-friendly palette for the color gradient
(with option = "magma" for the magma variant), using the
square root of the population (which is stored in the variable
POP_EST of the world object):
ggplot(data = world) +
geom_sf(aes(fill = pop_est)) +
coord_sf(expand = FALSE) +
scale_fill_viridis_c(option = "plasma", trans = "sqrt") # sqrt transform

coord_sf)The function coord_sf allows to deal with the coordinate
system, which includes both projection and
extent of the map. By default, the map will use the
coordinate system of the first layer that defines one (i.e. scanned in
the order provided), or if none, fall back on WGS84 (latitude/longitude,
the reference system used in GPS). Using the argument crs,
it is possible to override this setting, and project on the fly to any
projection. This can be achieved using any valid PROJ4 string (here, the
European-centric ETRS89 Lambert Azimuthal Equal-Area projection):
ggplot(data = world) +
geom_sf() +
scale_fill_viridis_c(option = "plasma", trans = "sqrt") +
coord_sf(crs = "+proj=laea +lat_0=52 +lon_0=10 +x_0=4321000 +y_0=3210000 +ellps=GRS80 +units=m +no_defs ")

Spatial Reference System Identifier (SRID) or an European Petroleum Survey Group (EPSG) code are available for the projection of interest, they can be used directly instead of the full PROJ4 string. The two following calls are equivalent for the ETRS89 Lambert Azimuthal Equal-Area projection, which is EPSG code 3035:
The extent of the map can also be set in coord_sf, in
practice allowing to “zoom” in the area of interest, provided by limits
on the x-axis (xlim), and on the y-axis
(ylim). Note that the limits are automatically expanded by
a fraction to ensure that data and axes don’t overlap; it can also be
turned off to exactly match the limits provided with
expand = FALSE:
ggplot(data = world) +
geom_sf(aes(fill = pop_est)) +
coord_sf(xlim = c(118, 128), ylim = c(17, 27), expand = FALSE) +
scale_fill_viridis_c(option = "plasma") # linear scale

ggspatial elementsSeveral packages are available to create a scale bar on a map
(e.g. prettymapr, vcd, ggsn, or
legendMap). Here the package ggspatial
provides easy-to-use functions. scale_bar that allows to
add simultaneously the north symbol and a scale bar into the
ggplot map. Five arguments need to be set manually:
lon, lat, distance_lon,
distance_lat, and distance_legend. The
location of the scale bar has to be specified in longitude/latitude in
the lon and lat arguments. The shaded distance
inside the scale bar is controlled by the distance_lon
argument while its width is determined by distance_lat.
Additionally, it is possible to change the font size for the legend of
the scale bar (argument legend_size, which defaults to 3).
The North arrow behind the “N” north symbol can also be adjusted for its
length (arrow_length), its distance to the scale
(arrow_distance), or the size the N north symbol itself
(arrow_north_size, which defaults to 6). Note that all
distances (distance_lon, distance_lat,
distance_legend, arrow_length,
arrow_distance) are set to "km" by default in
distance_unit; they can also be set to nautical miles with
“nm”, or miles with “mi”.
ggplot(data = world) +
geom_sf(aes(fill = pop_est)) +
coord_sf(xlim = c(118, 128), ylim = c(17, 27), expand = FALSE) +
scale_fill_viridis_c(option = "plasma") +
annotation_scale(location = "br", width_hint = 0.3) +
annotation_north_arrow(location = "br", which_north = "true",
pad_x = unit(0.5, "cm"), pad_y = unit(1, "cm"),
style = north_arrow_fancy_orienteering)

Note: if you plot a larger area, you may get a warning on the inaccuracy of the scale bar.
geom_text and
annotate)The world data set already contains country names and
the coordinates of the centroid of each country (among more
information). We can use this information to plot country names, using
world as a regular data.frame in
ggplot2. The function geom_text can be used to
add a layer of text to a map using geographic coordinates. The function
requires the data needed to enter the country names, which is the same
data as the world map. Again, we have a very flexible control to adjust
the text at will on many aspects:
The size (argument size);
The alignment, which is centered by default on the coordinates
provided. The text can be adjusted horizontally or vertically using the
arguments hjust and vjust, which can either be
a number between 0 (right/bottom) and 1 (top/left) or a character
(“left”, “middle”, “right”, “bottom”, “center”, “top”). The text can
also be offset horizontally or vertically with the argument
nudge_x and nudge_y;
The font of the text, for instance its color (argument
color) or the type of font
(fontface);
The overlap of labels, using the argument
check_overlap, which removes overlapping text.
Alternatively, when there is a lot of overlapping labels, the package
ggrepel provides a geom_text_repel function
that moves label around so that they do not overlap.
For the text labels, we are defining the centroid of the
countries with st_centroid, from the package
sf. Then we combined the coordinates with the centroid, in
the geometry of the spatial data frame. The package
sf is necessary for the command
st_centroid.
Additionally, the annotate function can be used to add a
single character string at a specific location, as demonstrated here to
add the “Gulf of Mexico”Pacific Ocean” and Ryukyu Archipelago
sf::sf_use_s2(FALSE) # FOR ERROR turn off the s2 processing
world_points<- st_centroid(world)
world_points <- cbind(world, st_coordinates(st_centroid(world$geometry)))
ggplot(data = world) +
geom_sf(aes(fill = pop_est)) +
geom_text(data= world_points,aes(x=X, y=Y, label=name),
color = "grey", fontface = "bold", check_overlap = FALSE) +
scale_fill_viridis_c(option = "plasma") +
annotate(geom = "text", x = 124, y = 21, label = "Pacific Ocean", fontface = "italic", color = "#0b3c8a", size = 5) +
annotate(geom = "text", x = 124.2, y = 24, label = "Ryukyu archipelago", fontface = "italic", color = "#d41919", size = 3) +
coord_sf(xlim = c(118, 128), ylim = c(17, 27), expand = FALSE)

Now to make the final touches, the theme of the map can be edited to
make it more appealing. We suggested the use of theme_bw
for a standard theme, but there are many other themes that can be
selected from (see for instance ?ggtheme in
ggplot2, or the package ggthemes which provide
several useful themes). Moreover, specific theme elements can be tweaked
to get to the final outcome:
Position of the legend: Although not used in this example, the
argument legend.position allows to automatically place the
legend at a specific location (e.g. “topright”, “bottomleft”,
etc.);
Grid lines (graticules) on the map: by using
panel.grid.major and panel.grid.minor, grid
lines can be adjusted. Here we set them to a gray color and dashed line
type to clearly distinguish them from country borders lines;
Map background: the argument panel.background can be
used to color the background, which is the ocean essentially, with a
light blue; Many more elements of a theme can be adjusted, which would
be too long to cover here. We refer the reader to the documentation for
the function theme.
ggplot(data = world) +
geom_sf(fill= 'antiquewhite') +
geom_text(data= world_points,aes(x=X, y=Y, label=name), color = 'darkblue', fontface = 'bold', check_overlap = FALSE) +
annotate(geom = 'text', x = -90, y = 26, label = 'Gulf of Mexico', fontface = 'italic', color = 'grey22', size = 6) +
annotation_scale(location = 'bl', width_hint = 0.5) +
annotation_north_arrow(location = 'bl', which_north = 'true', pad_x = unit(0.75, 'in'), pad_y = unit(0.5, 'in'), style = north_arrow_fancy_orienteering) +
coord_sf(xlim = c(-102.15, -74.12), ylim = c(7.65, 33.97), expand = FALSE) +
xlab('Longitude') + ylab('Latitude') +
ggtitle('Map of the Gulf of Mexico and the Caribbean Sea') +
theme(panel.grid.major = element_line(color = gray(.5), linetype = 'dashed', size = 0.5), panel.background = element_rect(fill = 'aliceblue'))

The final map is now ready, it is very easy to save it using ggsave. This function allows a graphic (typically the last plot displayed) to be saved in a variety of formats, including the most common PNG (raster bitmap) and PDF (vector graphics), with control over the size and resolution of the outcome. For instance here, we save a PDF version of the map, which keeps the best quality, and a PNG version of it for web purposes:
rgbifImport Global Biodiversity Information Facility (GBIF) data using rgbif
package and crete distribution map using the mapr
package.
gbif.res <- occ_search(scientificName = "Urocissa caerulea", limit = 1200)
map_ggplot(gbif.res) +
coord_sf(xlim = c(120, 123), ylim = c(21, 26), expand = FALSE)

Note the slight difference in overlap of the distribution data on this shape file.
Have a look here, to combine climate and species data.
marmapmarmap can query and plot NOAA’s bathymetry database
# query
TW.bathy <- getNOAA.bathy(lon1=118,lon2=124, lat1=21,lat2=26,resolution=1) # don't put too wide / resolution: 1
# define palette
blues <- colorRampPalette(c("darkblue", "cyan"))
greys <- colorRampPalette(c(grey(0.4),grey(0.99)))
# make the plot
plot.bathy(TW.bathy,
image=T,
deepest.isobath=c(-6000,-120,-30,0),
shallowest.isobath=c(-1000,-60,0,0),
step=c(2000,60,30,0),
lwd=c(0.3,1,1,2),
lty=c(1,1,1,1),
col=c("grey","black","black","black"),
drawlabels=c(T,T,T,F),
bpal = list(c(0,max(TW.bathy),greys(100)),c(min(TW.bathy),0,blues(100))),
land=T, xaxs="i"
)

Profiles can be extract using get.transect:
tw.profile <-get.transect(TW.bathy,x1=119.5,y1=23.75, x2=122,y2=23.75, dis=TRUE)
plotProfile(tw.profile)
#### Not Run: extract a profile Manually
#### manual.profile<-get.transect (TW.bathy, loc=T,dist=T)
#### plotProfile(manual.profile)
Leaflet is one of
the most popular open-source JavaScript libraries for interactive maps.
You can loose hours and its functionalities are amazing. The R package
leaflet makes it easy to integrate and control Leaflet maps
in R.
FRE <- paste(sep = "<br/>",
"<b><a href='https://www.dipintothereef.com/'>FRELAb TAIWAN</a></b>",
"Functional Reef Ecology Lab",
"Institute of Oceanography, NTU")
leaflet(taiwan) %>%
addPolygons(weight=0.5) %>%
addTiles(group="Kort") %>%
addPopups(121.53725, 25.021252, FRE, options = popupOptions(closeButton = FALSE))
⚠ Practice 3: Create and customize an interactive map
of your choice after exploring the functionality of the
leaflet package. Push both your .Rmd and
.html files into a public repository available from your Github
account.You will share with me be email [vianneydenis@g.ntu.edu.tw] the address (URL) of this
repository (such as https://github.com/vianneydenis/OCEAN5098.git)
before next Monday in order for me to check your work.
The title of your email should be `Practice 3 (your
name: your student no.). BE CREATIVE ;)
Text and figures are licensed under Creative Commons Attribution CC BY 4.0. The figures that have been reused from other sources don't fall under this license and can be recognized by a note in their caption: "Figure from ...".